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Abstract 

Due to the action of the intervening cosmic magnetic fields, ultra-high energy cosmic rays (UHECRs) can be deflected in such a 
PsJ way to create clustered energy-ordered filamentary structures in the arrival direction of these particles, the so-called multiplets. In 
jj^this work we propose a new method based on the spherical wavelet transform to identify multiplets in sky maps containing arrival 

directions of UHECRs. The method is illustrated in scenarios with a multiplet embedded in isotropic backgrounds with different 
^1 pumber of events. The efficiency of the algorithm is assessed through the calculation of Type I and II errors. 

Keywords: spherical wavelets, ultra-high energy cosmic rays, cosmic magnetic fields, multiplets 
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1. Introduction 

[ Cosmic rays were discovered more than one century ago. 
One remarkable feature of the cosmic ray spectrum is that it 
ppans more than ten orders of magnitude, up to hundreds of 
EeV (1 EeV =10^^^ eV). The spectrum roughly follows an in- 
verse power law, which means that the expected flux of particles 
at the highest energies is extremely low compared to the lower 
(energies. In fact, at energies of a few EeV, only one particle 
per square kilometer per year is expected. Particles with en- 
ergies > 1 EeV are referred to as Ultra-High Energy Cosmic 
Rays (UHECRs). Some experiments, such as the Pierre Auger 
Observatory and the Telescope Array, have been designed to in- 
crease the statistics of events in this energy range. Despite the 
^igher statistics achieved, which has shed light over many im- 
portant questions, some of them still remain unanswered, such 
as the origin, nature and mechanisms of acceleration of these 
particles. 

\ Due to the presence of galactic and extragalactic magnetic 
fields, charged cosmic rays are expected to be deflected. Hence, 
incoming directions, as measured by a detector, do not point 
back to the exact position of the source. The magnitude of 
the deflection depends on the strength of the intervening fields. 
In the case of charged particles the deflections are roughly in- 
versely proportional to the particle energy. Therefore, for co- 
herent fields, the different Larmor radii described by cosmic 
rays can create filamentary structures ordered by energy, the 
multiplets. This allows the reconstruction of the source posi- 
tion, and consequently enhances the possibility to do astronomy 
wifli UHECRs. 

In this paper we propose a new method of identifying mul- 
tiplets, based on the spherical wavelet transform. The paper is 



organized as follows: in section 2, we review the physics under- 
lying multiplets; in section 3 we present the wavelet transform 
from a pattern matching algorithm point of view, and give moti- 
vations for its use in cosmic ray physics; in section 4 we present 
a novel algorithm to identify filamentary structures in cosmic 
ray maps; in section 5 the method is applied to simulated data 
sets; in section 6, we present our results; in section 7 we make 
our final remarks. 

2. Cosmic Magnetic Fields and Multiplets 

Several results show that at ultra-high energies the cosmic 
ray spectrum might contain a component of atomic nuclei. Data 
from the High Resolutions Fly's Eye Experiment (HiRes) indi- 
cate that UHECRs are probably protons [1], whereas the results 
from the Pierre Auger Observatory indicate that the composi- 
tion tends to heavy nuclei at the highest energies' [2]. Despite 
this controversy, one can consider that UHECRs are predomi- 
nantly charged particles and, as so, can be deflected by mag- 
netic fields. 

The deflection expected for a UHECR of charge Z due to 
the regular component of the galactic magnetic field (GMF) is 
given approximately by [3]: 



E 



df B 
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EeV, 



(1) 



where E is the energy of the particle, B the magnetic field, and 
L the travel distance of the particle. Since \B\ ^ fiG and is 
coherent over lengths of ^ 10 kpc, typical deflections are ^ 
10° for protons of 10 EeV. 



Email addresses: mzimbres@if i . unicamp . br (M. Zimbres), 
rafael. alves. batistaSdesy . de (R. Alves Batista), 
kemp@ if i . unicamp . br (E. Kemp) 



Notice tliat tliis results strongly depends upon the hadronic interaction 
model taken into account. 
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The expected deflection due to the turbulent component is 



[4]: 



assume that the deflections induced by extragalactic magnetic 
fields are much smaller than the ones induced by the GMF. 
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(2) 3. Wavelets on the sphere 



where Lc is the coherence length of the field and Brms is the 
root mean square intensity of the magnetic field. 

According to equation (1) the deflection is inversely pro- 
portional to the energy of the particle. Therefore it is possible 
that energy ordered filamentary structures, the so-called multi- 
plets, can be detected in cosmic ray maps as shown in figure 1 . 
Another method to detect filamentary structures in the arrival 
directions distributions of UHECRs was proposed by Harari et 
al.[5]. 

An interesting property of the multiplets is that the posi- 
tion of the source can be reconstructed, allowing one to identify 
UHECRs sources. A method to reconstruct the source position 
of a multiplet was presented by Golup et al. [6] and was used 
by the Pierre Auger Collaboration to estimate the position of the 
sources for some possible multiplet candidates [4]. In the afore- 
mentioned work, no evidence for the presence of multiplets was 
found for energies above 20 EeV. 




Figure 1 : Illustration of a multiplet on the tangent plane of a 
sky map. The dots represent events from higher (red) to lower 
(blue) energies. The black dot corresponds to the position of 
the source. 

It is important to notice that for some models of the galactic 
magnetic fields such as the ones proposed by [7-9] cosmic ray 
multiplets can be formed, whereas in other models such as the 
one recently proposed by [10] they are very unlikely to occur, 
due to the strength of the field, especially the turbulent compo- 
nent. 

The role played by extragalactic magnetic fields in the de- 
flection of UHECRs is not fully understood. Simulations of 
the propagation of UHE particles in the large scale structure of 
the universe have been performed by several groups [11-14]. 
However, these results are contradictory and one cannot obtain 
a clear picture of the effects of the extragalactic magnetic field 
for the deflection of UHECRs. For the purposes of this work we 



In many branches of Physics, specially Astrophysics and 
Cosmology, wavelets have been successfully applied to solve 
various problems, particularly related to detection of signals. 
Wavelets on the plane have been widely used to denoise cos- 
mic microwave background (CMB) maps [15-17]. However, 
the problem of identifying anisotropics in the distribution of ar- 
rival directions of UHECRs has not been properly addressed, 
and only a few works [18-22] on this topic are available in the 
literature. 

Wavelets are commonly used in one and two dimensional 
data analysis, but in recent years the interest in data lying on 
the sphere has increased due to experiments such as the Cos- 
mic Background Explorer (CoBE), the Wilkinson Microwave 
Anisotropy Probe (WMAP), the Planck Satellite and the Pierre 
Auger Observatory, which make use of these kind of data. The 
need to process these sorts of data sets pushed the interest on 
new techniques under development, being wavelet analysis on 
the sphere one example. Wavelets on the sphere have been 
successfully applied to studies of the cosmic microwave back- 
ground [23-26], and also dark energy detection [27, 28]. 

Wavelets can be specially useful for cosmic ray data analy- 
sis, where we usually have to deal with a non uniform exposure 
to search for local structures in the sky, with a defined position, 
such as point sources, and possibly an orientation, such as mul- 
tiplets. An event by event analysis is not viable, and the analysis 
in harmonic space would be even harder since all local proper- 
ties are lost. Spherical wavelets come up as a good alternative 
to address these problems. Other attractive features of wavelet 
analysis are: 

• any function can be exactly represented by its wavelet 
coefficients; 

• local features of the signal can be enhanced in wavelet 
domain, meaning that the number of coefficients needed 
to represent a given signal is reduced; 

• it provides scale decomposition, making it possible to 
identify structures with different angular sizes and focus 
on the part of interest in the signal; 

• it does not rely on any tangent plane approximation, in 
the case of wavelets on the sphere. 

Data representation in wavelet domain can be thought as 
something between pixel and harmonic representation. Some- 
times it is very convenient to decompose the data to enhance 
properties that are not clear in harmonic or pixel domain. Wavelets 
can be interpreted as local analysis functions which can be ro- 
tated and/or dilated, to obtain information regarding the signal 
morphology. 
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3.1. Pattern matching on the sphere 

In this subsection we show that the problem of finding a 
multiplet, or any other pattern defined on the sphere, can be 
treated by the fast rotational matching algorithm [29]. 

Let f{9,(p) and h{9,(p) be two functions defined on the 
sphere. Assume that ft, is a rotated version of /, such that 
/ = A(a, /3, ^)h, with (a, /?, 7) being Euler angles, and A de- 
noting the rotation operator in S'0(3). 

If we know that a rotated version of the pattern / is present 
in h, we can find its latitude, longitude and orientation on the 
sphere by correlating all rotated versions of / with h, and se- 
lecting the rotation which maximizes the correlation 

c = j h{e,(p)Af{e,(p)d{cose)dip. (3) 

We can parametrize the rotations in terms of Euler angles. In 
this case the correlation function can be denoted by C = C{a,/3,'^). 
In other words, we want to find the angles a, (3 and 7 for which 
C is maximum. 

The straightforward evaluation of C is very time-consuming 
and not affordable depending on the precision we are using to 
describe our signal. Denoting the band-limit of the signal by 
B, the complexity of equation (3) is 0{B^). This complex- 
ity can be reduced to 0{B^) by calculating the correlation in 
the harmonic domain. So, we can write the spherical harmonic 
expansions of / and h as 

B-l I 

f{e,ip)^Y. E «'™^™(^'^) (4) 

B-l I 

h{e,^)^J2 E biruYiie,^), (5) 

with YJ^ being spherical harmonics. Using these expansions we 
can show (see Ref. [30] for more details) that the correlation 
function (3) can be cast as 

C = E E 0^^lnD''mn{ot^Pn) (6) 
m,n— — ^ 
B-l I 

= J2 J2 «„e^™"dL(/5)e^"^ (7) 

/— m.n— — l 

where -Dmn denotes the Wigner-D functions and dl^„ the small 
Wigner-d functions. This algorithm is called Fast Rotational 
Matching, abbreviated to FRM. It uses the Fast Fourier Trans- 
form (EFT) on both a and 7 to calculate equation 7. Depending 
on the tessellation chosen, the EFT can be extended to the f3 
coordinate [29]. 

This approach to find a multiplet may not be an option be- 
cause we do not know exactly the pattern f{6, ip) that describes 
the multiplet. Moreover, for angular resolutions corresponding 
to a band limit higher than B = 128, too much memory would 
be required. For B = 256, for example, approximately 4.5 GB 
of RAM would be necessary in our implementation. 

In the next subsection we show that the difficulties men- 
tioned above can be overcome by using a special family of di- 
rectional wavelets instead of the pattern itself. 
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Figure 2: The figure shows the harmonic support in the fre- 
quency domain of the family of wavelets used, for J ~ 8. 

Table 1: First column: scale j. Second column: wavelet sup- 
port. Third column: angular sizes in degrees for which the 
wavelet is sensitive. Fourth column: maximum precision on 
the angular variables in degrees. 



j Support Angular size (°) Precision (°) 






(256,128) 


(0.7,1.4) 


0.7 


1 


(256,64) 


(0.7,2.8) 


0.7 


2 


(128,32) 


(1.4,5.6) 


1.4 


3 


(64, 16) 


(2.8,11.3) 


2.8 


4 


(32,8) 


(5.6,22.5) 


5.6 


5 


(16,4) 


(11.3,45.0) 


11.3 


6 


(8,2) 


(22.5,90.0) 


22.5 


7 


(4,1) 


(45.0,180.0) 


45.0 


8 


(2,1) 


(90.0,180.0) 


90.0 



3.2. Pattern matching with directional wavelets 

When the function / from equation 4 is a wavelet, and h is 
the signal of interest, equation 3 is the definition of the forward 
spherical wavelet transform, and the correlation function can be 
interpreted as the wavelet representation of the signal h. 

The wavelets used in this work are defined in the harmonic 
domain^. Their harmonic representation has the special prop- 
erty of allowing them to be split into a kernel and a directional 
part 

blm = k{l)Slm, (8) 

where the kernel k{l) is responsible for dilations and Sim is 
responsible for the directional properties of the wavelets. This 
spht ensures that dilations do not affect directional properties, 
so that these two parts can be treated independently. 



^For the explicit definition and derivation of the mother-wavelet equation, 
refer to Ref. [31]. 
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The kernel k{l) at each scale j has zero values at frequen- 
cies outside the range (2''^^^-', 2'^+^^-'), which is usually re- 
ferred to as the wavelet support. The support is related to the 
frequencies to which the wavelets are sensitive. For a graphical 
representation see figure 2. It is interesting to notice in this fig- 
ure that an adequate choice of j can suppress some values of I, 
making this method powerful even when we consider the expo- 
sure of a detector, which is usually associated to low values of 

In our implementation, for a given band limit B, the num- 
ber of scales in the wavelet analysis is given by J = log2 B. 
Since we have used B ~ 256, then J = 8. For a physical in- 
terpretation, however, it is simpler to think in terms of angular 
sizes. For that, we can convert the frequency range into an an- 
gular size using the formula 360°/2/, where I is the frequency 
we are interested in. In table 1 the angular size to which each 
scale is sensitive is shown. 

The sensitivity of the wavelet on finding the angle 7 can be 
controlled by imposing a band limit on Sim- Denoting this band 
Umit by N, then Sim = for all m > N. If the directional fea- 
tures of the signal are not relevant for the analysis, low values 
of N, such as = 1, can be used. The maximum value of N 
at scale j is given by = 2'^^^+^. The precision of the orien- 
tation for a given N is given by A7 = 180/A^. In this work, we 
have used N = 127, shown in figure 3, which gives a precision 
A7 = 1.42 degrees. This value is a good compromise between 
computational resources and the required precision. 




l -in I- 

-30 -20 -10 10 20 30 -30 -20 -10 10 20 30 



(a) TV = 1 (b) N = 127 

Figure 3: Comparison of wavelet with parameters ( J, j, N) = 
(8, 2, 1) and ( J, j, N) = (8, 2, 127). Both axes are in units of 
degrees. The wavelet shown in figure (b) is the one used in the 
analysis. 



4. Looking for multiplets 

We have developed an algorithm that uses spherical wavelets 
to identify and locate filaments in maps containing arrival direc- 
tions of UHECRs. The algorithm is described below. 

1. From a set of events, calculate the function h{d, ip) that 
represents the signal. This involves choosing a tessel- 
lation for the sphere, and counting the total number of 
events contained in each pixel in the sky. 

2. Calculate the Fourier expansion of h{d, ip), resulting in 
the coefficients aim (see equation 4). 



3. Choose the appropriate wavelet. For the family of wavelets 
used, this involves choosing three parameters, the maxi- 
mum scale J, the scale at which the analysis is performed 
j and the azimuthal band limit, controlled by the parame- 
ter N. Ideally, the angular size of the wavelet will match 
the size of the multiplet. This choice relies on 1 . 

4. Calculate C{a, (5, 7) from equation (7). 

5. Select all (a, /3,7) such that |C(a, /3,7)| > Co, where 
Co, is threshold value. 

6. Select the events at each location (a,/3,7) (see section 
4.2). This will result in m groups of events. 

7. Discard all groups of events for which n < jiq, where n 
is the number of events in the group, and uq is a threshold 
value. 

8. For each group of events, calculate the correlation c of 
the graph S x 1/E (see equation 17). 

9. Accept as a multiplet candidate all groups for which |c| > 
Co, where cq is a threshold value, (see section 4.2) 

4.1. Establishing thresholds 

To carry out all steps of this algorithm one needs to establish 
three threshold values: the wavelet threshold Co used in step 
5, the minimum number of events no, used in step 7, and the 
minimum correlation cq, used in step 8. 

To calculate the thresholds, we have used simulated sky 
maps containing arrival directions of events isotropically dis- 
tributed. No magnetic fields are considered in this case, so that 
if a multiplet is identified by our method, it certainly happened 
by chance. This allows us to establish the threshold values Co, 
Co, and uq by using algorithm previously presented in section 
4. 

To estimate Co, we start with M isotropic realizations fol- 
lowing the same assumptions on the injection spectrum and ex- 
posure. Ci is the largest wavelet coefficient obtained for each 
realization i, according to step 5 of the algorithm. 

By selecting only one value of C for each realization, we 
have a single correlation coefficient c, for each of the simulated 
isotropic skies, calculated at step 8 of the algorithm. 

From M wavelet coefficients Ci and correlations Ci, we 
choose 

Co = C + rccrc and co = c + rc<Jc, (9) 

where C and c denote the average values over M realizations, 
and (7(7 and CTc their respective standard deviations. The num- 
bers rc and rc are chosen according to the error the user is 
wilUng to accept. In this work we set = 1 and rc ~ —1. 

The threshold value for the number of events in the multi- 
plet, no, is chosen depending on the specific problem. In this 
analysis, for example, we have chosen uq = 10, following Ref. 
[4]. 

4.2. Selecting events 

To select events around the position where a wavelet has a 
high magnitude we need the Euler angles a and /3, to provide a 
location in the sky, and the angle 7, to provide an orientation. 
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A multiplet can, in principle, describe an arbitrary curve on the 
surface of the celestial sphere, depending on the intervening 
magnetic fields. To be general for all possible shapes we have 
used small segments as shown on figure 4. Each segment should 
be small enough so that the curve described by the events is 
approximately a straight line inside the segment. 

A small deflection angle (< 10°) can be written as 6 = 
As /r, which in spherical coordinates takes the form 

As^ = r'^{{A/3f + sin^ l3{Aaf), (10) 

with As being the linear distance over the surface of a sphere 
of radius r. 

Now we note that the unit vectors a and /3, which point in 
the direction where a and /3 vary, are orthogonal and span the 
plane tangent to the sphere at (a, f3). Using equation 10 we can 
write a displacement vector, tangent to that point as 

d = r{A(3)$ + rsm(3{Aa)a. (11) 

It is convenient to write d in terms of L, aligned with the wavelet, 
and W, which is perpendicular to it. These vectors are related 
to a and /3 by a rotation of an angle 7 in the tangent plane, as 
shown below: 

a = cos 7L — sin 714^ (12) 
/3 = sin7i7 — cos7W^. (13) 

We can now cast 1 1 as 

d:^diL + d^W, (14) 

where 

dj^ = r(A/?cos7 + sin/3Aasin7) (15) 
d^ ~ r(— A/3sin7 + sin/3AQ;cos7). (16) 

With equations 15 and 16, we can select events in the sky at 
the position of maximum wavelet coefficient. The only infor- 
mation required are the Euler angles a, (3, and 7, since r = 1 
for the celestial sphere. 

5. Analysis 

To simulate both the background and the events belong- 
ing to the multiplet, we have used the CRPropa software [32]. 
Protons were injected following a typical E^^-"^ spectrum, for 
15 < i? < 40 EeV. The effective detector size is a sphere of 
radius 1 Mpc, and the source was assumed to be in our local 
universe, at a distance of 30 Mpc from the detector For the 
isotropic data sets, the events were simulated in random posi- 
tions of the sky. 

We have divided our analysis in two parts. First we analyze 
an isotropic distribution of events to establish the thresholds, as 
explained in section 4. 1 . In the second part of the analysis we 
used a simulated multiplet composed by 10 events, embedded 
in the same isotropic datasets previously used. The multiplet 
formation was induced by an unrealistic uniform magnetic field 
of 1 nG oriented in the z direction. Nevertheless, the results are 




Figure 4: Illustration of a multiplet along a segment. The dots 
represent events from higher (red) to lower (blue) energies . The 
angles (a, /3, 7) are the Euler angles associated with the multi- 
plet. 



totally independent of the method used to generate the multi- 
plets, since it depends only on the position of the events, and 
on the energy- deflection correlation. The number of isotropic 
events range between 100 and 1000 events, in steps of 100, and 
a total of 1000 realizations for each of these values. The mul- 
tiplet embedded in one of the isotropic skies containing 1000 
isotropically distributed events can be seen in figure 5. 




20 25 30 35 

Figure 5: Simulated multiplet containing 10 events embedded 
in a background of 1000 events. The color scale corresponds to 
the energy of the events in units of EeV. The zoom shows more 
closely the events of the simulated multiplet. 
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The simulated multiplet has a correlation coefficient c = 
0.99 when no background is present, c corresponds to the Pear- 
son's coefficient of the 6 x E^^ graph, given by: 

c = '"i (17) 

The scale at which we performed the analysis was chosen 
based on table 1. At scales j = 2 and 3 the angular size of 
the wavelet is close to what we would expect for a multiplet 
similar to the simulated one, i.e., 10° x 2°. Since at scale j = 
3 the maximum precision of the angular variable is 2.8°, we 
have used j = 2 for the analysis. The size of the segments are 
10° X 2°, to match the dimensions of the wavelet. 

6. Results and Discussion 

In this section we present the results of the proposed algo- 
rithm, when applied to the simulations. 

Since the magnitude of the wavelet coefficients is used to 
select directions of interest in sky maps, i.e., Euler angles, es- 
tablishing the wavelet coefficient threshold is a crucial part of 
the analysis. The wavelets should have enough sensitivity to 
distinguish the cases of interest (isotropic data sets with an em- 
bedded multiplet and without it). Figure 6 shows the aver- 
age values of the magnitude of the largest wavelet coefficients 
found in each realization, for the two data sets. The error bars 
are the standard deviation of the 1000 realizations of isotropic 
skies. In this particular case, for a multiplet composed of 10 
events, the error bars start overlapping at around 1000 events, 
which corresponds to a fraction of events from the multiplet 
with respect to the background of w 10^^. 

The correlation coefficient is the main observable to prop- 
erly characterize a multiplet. In this analysis we have used 
no = 10, which means that a candidate will only be considered 
a multiplet if it is composed by, at least, 10 events. We also have 
to determine the value of co, which is the value of the thresh- 
old correlation coefficient. It is expected that the greater the 
number of isotropically distributed events, which corresponds 
to the background, the smaller the correlation coefficient. This 
is corroborated by figure 7, which shows that the correlation 
coefficient monotonically decreases for larger number of back- 
ground events, indicating a contamination of the multiplet by 
the background. This allows us to establish a threshold value 
for the correlation coefficient. For 1000 isotropic events the 
lowest value of cq considered is cq = 0.4, according to equa- 
tion 9. 

It is important to address the questions related to the proba- 
bility of the algorithm finding a multiplet with c > co in a given 
data set, when no multiplets are known to be present (Type II 
error), and the probability of a multiplet to have c < cq if it is 
known to be present in the data (Type I error). In the isotropic 
simulations, no multiplets with c > 0.4 were identified for a 
fraction of events inferior to 10^^, implying that the maximum 
type II error of the method is lO^'^. Figure 8 shows the Type 



Magnitude of wavelet coefficients |C| 
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Figure 6: Magnitude of the largest wavelet coefficients in all 
the 1000 realizations, as a function of the number of events in 
the background. Blue circles correspond to the purely isotropic 
case, and red squares to the multiplet embedded in the isotropic 
background. The error bars are the standard deviation of the 
wavelet coefficients for all the realizations. 



Correlation comparison 




200 400 600 800 1,000 
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Figure 7: Correlation coefficients of the 1000 realizations, as 
a function of the number of events in the background. Blue 
circles correspond to the purely isotropic case, and red squares 
to the multiplet embedded in the isotropic background. The 
error bars are the standard deviation of the wavelet coefficients 
for all the realizations. 



I error introduced by our method. It is worth mentioning that 
for the most conservative case, which has the lowest fraction of 
number of events of the multiplet with respect to the number of 
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events of the background, the Type I error introduced by using 
a threshold value cq = 0.4 is approximately 25%, whereas for 
the least conservative case it is around 2%. For a high threshold 
value such as cq = 0.9, this error is of the order of 30% for 
the most conservative case. It is important to stress that the es- 
timated Type I and II errors depend on the confidence interval 
adopted, set by the parameters rc and Tc. 
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Figure 8: Probability of having c < cq (type I error) when the 
multiplet is known to be present in the data set. Each curve 
corresponds to a different fraction of events of the multiplet 
with respect to the isotropic background. 



7. Conclusion 

We have developed a wavelet based analysis method to iden- 
tify multiplets in maps containing arrival directions of UHE- 
CRs. We have illustrated the method by applying it to a hypo- 
thetical scenario with a uniform magnetic field of 1 nG, consid- 
ering several fractions of events from the multiplet with respect 
to the background, down to a fraction of 10^^. All the param- 
eters used in this study are rather arbitrary. These choices were 
made to illustrate the method, and not to evaluate its overall 
performance. 

The observables used to accept or reject the multiplet candi- 
date were the magnitude C of the wavelet coefficient at the mul- 
tiplet location (a, /?, 7), the correlation c of the deflection ver- 
sus inverse of the energy graph and the total number of events 
in the multiplet uq. 

The efficiency of the method was assessed by comparing the 
analysis when applied to isotropically distributed events with 
and without an embedded multiplet. The probability of wrongly 
accepting a candidate multiplet with a correlation coefficient c 



above the threshold value cq > 0.4 is below 10"'^. The prob- 
ability of wrongly rejecting a candidate multiplet, when it is 
known to be present in the data, is in the most conservative case, 
approximately 30%. This value goes down to 3% by decreasing 
the number of isotropic events in the data set. For the adopted 
correlation coefficient threshold, these values are, respectively, 
25% and 2%. 

Even though we have assumed a uniform detector exposure 
for the analysis previously presented, the results also hold for a 
non uniform exposure, since the ideal parameters for multiplet 
search in the wavelet analysis suppress low frequency modes, 
that could be related to the coverage of the cosmic ray detec- 
tor The only caveat is that, in this case, the comparison data 
sets should follow an isotropic distribution modulated by the 
exposure of the detector 

The method focused on the search of multiplets, but it can 
be adapted for generic application to related problems in astro- 
physics, particularly the ones involving the search of filamen- 
tary structures in sky maps. 

All results showed in this paper can be easily reproduced 
with the software SWAT (Spherical Wavelet Analysis Tool), de- 
veloped by the authors of this paper. In case of interest in the 
code, please contact us. 
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